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Properties of the one-dimensional totally asymmetric simple exclusion process (TASEP), and 
their connection with the dynamical scaling of moving interfaces described by a Kardar-Parisi- 
Zhang (KPZ) equation are investigated. With periodic boundary conditions, scaling of interface 
widths (the latter defined via a discrete occupation-number-to-height mapping), gives the exponents 
a = 0.500(5), 2 = 1.52(3), f3 = 0.33(1). With open boundaries, results are as follows: (i) in the 
maximal-current phase, the exponents are the same as for the periodic case, and in agreement with 
recent Bethe ansatz results; (ii) in the low-density phase, curve collapse can be found to a rather 
good extent, with a = 0.497(3), z = 1.20(5), j3 = 0.41(2), which is apparently at variance with 
the Bethe ansatz prediction z = 0; (iii) on the coexistence line between low- and high- density 
phases, a = 0.99(1), z = 2.10(5), (i = 0.47(2), in relatively good agreement with the Bethe ansatz 
prediction z = 2. From a mean-field continuum formulation, a characteristic relaxation time, related 
to kinematic- wave propagation and having an effective exponent z = 1, is shown to be the limiting 
slow process for the low density phase, which accounts for the above-mentioned discrepancy with 
Bethe ansatz results. For TASEP with quenched bond disorder, interface width scaling gives a — 
1.05(5), 2 = 1.7(1), B — 0.62(7). From a direct analytic approach to steady-state properties of 
TASEP with quenched disorder, closed-form expressions for the piecewise shape of averaged density 
profiles are given, as well as rather restrictive bounds on currents. All these are substantiated in 
numerical simulations. 

PACS numbers: 05.40.-a,02.50.-r,05.70.Fh 



I. INTRODUCTION 

In the absence, so far, of a general theory describ- 
ing non-equilibrium processes (even in steady state), it 
is worthwhile studying simple models whose qualitative 
features, it is hoped, will hold for a broad class of sys- 
tems. In this paper we deal with properties of the one- 
dimensional totally asymmetric simple exclusion process 
(TASEP) fjj, and their connection with the dynamical 
scaling of moving interfaces described by a Kardar-Parisi- 
Zhang (KPZ) equation 0,0,11. The TASEP is a biased 
diffusion process for particles with hard-core repulsion 
(excluded volume) 0,H, @|. 

Our main purpose is twofold: first, to probe the re- 
lationship between TASEP behavior and KPZ interface 
evolution under several distinct constraints, to be de- 
scribed below; and, focusing especifically on systems with 
quenched disorder, to provide an account of the effects 
of frozen randomness on the density profiles and currents 
in TASEP. 

Upon implementing specific features on the particle 
system, such as various boundary conditions, assorted 
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particle densities, currents, and/or injection/ejection 
rates, as well as quenched inhomogeneities, we measure 
the consequent changes to properties of the interface 
problem, taking the latter to be related to the former 
by the connection which we now sketch. 

The 1 + 1 dimensional TASEP is the fundamental dis- 
crete model for flow with exclusion. Here the particle 
number ni at lattice site £ can be or 1, and the for- 
ward hopping of particles is only to an empty adjacent 
site. Taking the stochastic attempt rate p = 1 (see 
later) the current across the bond from I to £ + 1 is 
thus Jtij+i = ng(l — n/ + i). This system maps exactly 
to an interface growth model Q in D = 1 + 1, having 
integer values of height variables h^p, on a new lattice 
such that lies midway between sites £, I + 1 of the 
TASEP lattice. The hi are constrained by the relation 
h>M+i) ~ him = 1 - 2rif (= ±1). When not concerned 
with detailed associations of sites between the models we 
will omit the (£) in i{£); after continuum limits, {£) and 
i{£) become the same. 

There is clearly a precise one-to-one correspondence 
between these two discrete models and their stochasticity. 
Various associated continuum models are not so clearly 
related, due to the possibility of different continuum lim- 
its. A naive continuum limit in the TASEP gives 0, 0] 
the noisy Burgers turbulence equation 0, [n], EH for one- 
dimensional particle flow: this equation is the continu- 
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ity equation resulting from a continuum version of the 
TASEP bond current J, which takes the form 

'~5£-H) > + «*" + i' (l) 

where p = (n) (i.e., the "mean field" version of n), and 
t) is the (uncorrelated) noise here used to represent 
all the effects of stochasticity: (((x,t) ((x' ,t')) = 5(x — 
x')S(t-t'). 

Using the continuum "noiseless" version of the 
height/occupation given above for the discrete models, 



dh 
dx 



= 1 - 2p , 



(2) 



it is seen that the Burgers equation is related by a spatial 
derivative to the KPZ equation [2] for the evolution of 
the height h(x,t) of an elastic interface above a fixed 
reference level: 



dh 



dx 2 



Oh 

dx 



+ C(x,t) . 



(3) 



Depending on boundary conditions imposed at the 
extremities, one can have different regimes for the 
TASEP whose features will be recalled below where 

pertinent. There should be corresponding regimes also 
for KPZ, as pointed out previously [l2l.[l3|. 

In Section [ill we investigate the TASEP with periodic 
boundary conditions, and the associated interface prob- 
lem. By examination of interface width scaling, we es- 
timate the respective critical exponents. We also study 
the behavior of asymptotic interface widths against par- 
ticle density in the TASEP, as well as the main features 
of interface slope distributions and their connections to 
TASEP properties. In Sec. IIII1 we examine the interface 
width evolution corresponding to open-boundary TASEP 
systems in the following phases: (i) maximal-current, (ii) 
low-density, and (iii) on the coexistence line. In the lat- 
ter case, we also calculate density profiles in the particle 
system. In Sec. |IV] we turn to quenched bond disorder 
in the TASEP with PBC, providing a scaling analysis of 
the associated interface widths, as well as results for in- 
terface slope distributions. A direct analytic approach 
to steady-state properties of the TASEP with quenched 
disorder is given, focusing on averaged profiles densities 
and system-wide currents. Finally, in Sec. [VJ concluding 
remarks are made. 



II. PERIODIC BOUNDARY CONDITIONS 

We start by imposing periodic boundary conditions 
(PBC) for the TASEP at the ends of the chain, thus 
the total number of particles is fixed. Henceforth, the 
position-averaged particle density (p) will be denoted 
simply by p. Several steady-state properties are known 



exactly in this case as the configuration weights are 
factorizable at stationarity. To make contact with KPZ 
interface properties, we consider the width w(L,t) of an 
evolving interface of transverse size L in 1 + 1 dimensions 
in the discrete height model outlined in Section |H The 
average interface slope is 1 — 2p [ see Eq. ], and this 
average tilt must be taken into account. This is done by 
defining 

L 

[w(L : t)} 2 = L- 1 J2(h t (t)-h 1 (t)-(l-2p)(t-l)f , 

i=l 

(4) 

where only fluctuations around the baseline trend are 
considered. Initially, w grows with time as t°, until 
a limiting, L-dependent width ~ L a is asymptotically 
reached. With z = a/ (3, one expects from scaling 0,0]: 



w(L,t) = L a f 



L* ' 



where 



/(«) 



J vP u < 1 
1 const. 



(5) 



(6) 



For the D = 1 + 1 KPZ model, one has the exact val- 
ues 0: a = 1/2, (3 = 1/3, z = 3/2. 

We have simulated the TASEP on lattices with L — 
325, 650, and 1300 sites with PBC. For specified den- 
sities p, we would start from a particle configuration 
as uniform as possible, in order to minimize the asso- 
ciated interface width. A time step is defined as a set 
of L sequential update attempts, each of these according 
to the following rules: (1) select a site at random, and 
(2) if the chosen site is occupied and its neighbor to the 
right is empty, move the particle. Thus, the stochastic 
character resides exclusively in the site selection process. 
At the end of each time step we measured the width 
of the corresponding interface configuration. We took 
typically N s — 10 4 independent runs, averaging the re- 
spective results for each t. The evolution of interface 
widths for p = 1/2 and assorted lattice sizes is shown 
in Fig. [TJ the goodness of their scaling with the exact 
KPZ indices is typical of what is attained over the full 
density range. By examining the variation of data col- 
lapse quality against changes in the fitting exponents, we 
estimate a = 0.500(5), z = 1.52(3) for p = 1/2. Direct 
measurement of the short-time exponent (3 is less accu- 
rate; for example, fitting the L = 1300 data of Fig. Q] for 
10 2 < t < 10 3 to a single power law gives j3 = 0.31(1). 

The limiting (asymptotic) interface width, wy lm obeys 
the p <-» 1 — p particle-hole symmetry of the TASEP, 
with a maximum at p = 1/2. Fig. [2] exhibits the be- 
havior of w\i m for p < 1/2, and lattice size L = 325. 
The leftmost data point corresponds to eight particles, 
i.e. p = 0.0246. For even lower densities, discrete-lattice 
effects become more prominent, and the time needed to 
attain asymptotic behavior increases significantly. 



3 




Figure 1: (Color online) Double-logarithmic plot of interface 
width against time, corresponding to TASEP with periodic 
boundary conditions, p = 1/2, lattice sizes L as in key to 
symbols. Full line corresponds to w ~ t 1 ' 3 . Inset: scaling 
plot of data on main diagram, using a = 1/2, z = 3/2; L* = 
L/325. 



From Eq. (@]), using /ij — hi-i = 1 — 2ra J ;_i, the squared 
width is 



Hi,t)f = -^J( 1 -i) /) -^ nt 



(7) 



m — 1 



For the steady state with PBC, the configurational 
weights are factorizable [H]: each n m is indepen- 
dently distributed, taking values (0,1) with probabili- 
ties (p, 1 — p) . Consequently, in this case the average of 

(« ~ !)P-Emi n " 

state rms width is 

2\1 V 2 



is (i — 1 ) (p — p 2 ) , and the steady 
c(L-l)Va[p(i_ p )]i/a , (8) 



where c is a numerical constant. This form, consistent 
with having the KPZ scaling exponent a = 1/2, is com- 
pared with the simulation results for averaged limiting 
widths in Fig. [21 

According to Eq. ^ , particle density fluctuations can 
be investigated via the probability distribution functions 
(PDF) for slopes of the associated interface problem. In 
doing so within the context of a discrete-lattice model, 
one must take recourse to a coarse-grained description. 
We have experimented by taking interface segments with 
a varying number to of bonds, and calculating the aver- 
age slopes between the respective endpoints. For to > 10 
(the smallest practicable limit, such that the discrete- 
ness of allowed slope values still permits one to speak of 



Figure 2: (Color online) Limiting interface width wjim against 
particle density p for lattice with L = 325, PBC. Squares are 
numerical results, from N B — 10 4 samples each. Error bars 
are smaller than symbol sizes. The full line is the analytic 
expression, Eq. |[8|), with a = 9.97 (see text). 



a relatively smooth PDF), the curves are very close to 
Gaussians, and their width (height at peak) varies with 
to as to -1 / 2 (to 1 / 2 ). 

This can be understood by recalling the specific form 
of the_( factorizable) particle configuration weigths in this 
case 11. Taking p as the overall density, the probability 
of occurrence of a configuration with average slope s = 
1 — 2x (i.e., average local density x £ [0, 1]), on a lattice 
section with to sites is 



P{x) ~ cz x P mx 



i(l-x) 



(9) 



Standard treatment of Eq. ((9|) shows that, close to the 
maximum at x = p, the curve shape is indeed Gaussian, 
with a width proportional to to -1 / 2 . 

Having thus established the nature of the in- 
dependence of slope PDFs, we have used m = 60 in our 
calculations, which gives a convenient, almost continu- 
ous, spectrum of allowed slopes. Results for the station- 
ary state are depicted in Fig.[3l One sees that the p = 0.5 
PDF is fitted by a Gaussian down to some four orders of 
magnitude below the peak, while for p — 0.2, departures 
from a Gaussian profile are noticeable already at PDF 
values ~ 10 -2 times those at the maximum. 

We have followed the evolution of interface slope PDFs 
during the transient regime (starting from a particle dis- 
tribution as uniform as possible), for p — 1/2. We as- 
certained that, already from early times, their shapes 
are very well approximated by Gaussians. Fig. |4] shows 
the (root-mean-square) widths of Gaussian fits to the 
PDFs against time. It is noteworthy that, contrary to 
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the behavior of interface widths shown in Fig. [H here 
one does not find any significant dependence on lattice 
size. It thus appears that, even during the transient, 
the range of density fluctuations for p = 1/2 is shorter 
than the smallest lattice size considered here, L = 325. 
The solid curve in Fig. [4] is a stretched exponential, 
F(t) — a — b exp { — (t/to) s }, for which the best-fitting 
parameter values are to = 380(20), 5 = 0.51(2). 

In summary, we have shown that our methods, namely 
inferring scaling properties of TASEP via those of the 
associated interface problem, do give rather accurate re- 
sults (where comparison is possible, i.e., for the scal- 
ing exponents a, /3, and z) in the simplest case of a 
purely stochastic system with PBC. Our results so far 
are in agreement with the so-called KPZ conjecture 
sketched in Eqs. J!])-© above. 



III. OPEN BOUNDARY CONDITIONS 



Figure 3: (Color online) Slope PDFs in stationary state, lat- 
tice size L = 325 with PBC, for p = 1/2 and 0.2. Full lines 
are Gaussian fits to data. 



0.12 
0.1 
0.08 
0.06 
0.04 
0.02 



-III 


i 


1 1 










: <yo.07 




1 1 Mj. 




~Jf 0.06 






L= — 


4 0.05 






325 □ : 


f 0.04 






650 a j 








1300 o I 


6 0.03 


r i 


, I 1= 




I ° 

i i i i i 


200 \ 

i i 1 i i i 



2000 4000 

t 



Figure 4: (Color online) Width a of slope distributions in 
transient regime against time, for p = 1/2 and assorted lat- 
tice sizes L with PBC. Initial particle configuration: alternate 
empty and filled sites, for minimum interface width. Full line 
is stretched exponential fit to data (see text). Inset shows 
details of main figure, close to the vertical axis. 



A. Numerical results 

With open boundary conditions, the following addi- 
tional quantities are introduced: the injection rate aj at 
the left end, and the ejection rate (3e at the right one 
(both defined as fractions of the internal hopping rate) . 
The number of particles is no longer constant, although 
at stationarity it fluctuates around a well-defined aver- 
age. Therefore, in order to consider the associated in- 
terface widths, one needs to subtract the instantaneous 
average slope 1 — 2p(t), in the manner of Eq. (@]), at each 
time step. 

Many stationary properties are known for this case [l|, 
S S [3; EH, [IB|, including the phase diagram in 
(a/, Pe) space. With open boundary conditions one must 
be aware that, even at stationarity, ensemble- averaged 
quantities such as densities will locally deviate from 
their bulk values, within "healing" distances from the 
chain extremities which depend on the boundary (injec- 
tion/ejection) rates. 

Numerical density-matrix renormalization group 
(DMRG) techniques have been applied to estimate the 
dynamic exponent z for several locations on the phase 
diagram [T3|. More recently [3], a Bethe-ansatz solution 
has been provided, giving exact (analytic) predictions 
for the value of z everywhere on the phase diagram. 

We started by investigating the maximal-current (MC) 
phase (J = 1/4) at a 7 > 1/2, (i E > 1/2. The Bethe 
ansatz solution [13| predicts the KPZ value z = 3/2 
there, concurrently with DMRG results We took 

&i = Pe = 3/4, starting from an initial configuration 
with p = 1/2 (in the present case, this is the average fi- 
nal density as well), and alternate empty and filled sites, 
for minimum interface width. The evolution of interface 
widths against time was very similar to the PBC case, 
as shown in Fig. O and the KPZ exponents a = 1/2, 
/3 = 1/3, z = 3/2 were extracted within error bars of 
the same order as those for PBC. We also checked the 
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Figure 5: (Color online) Double-logarithmic plot of inter- 
face width against time, corresponding to TASEP with open 
boundary conditions inside the maximal current (MC) phase, 
cei = Pe = 3/4. Initial density p = 1/2, lattice sizes L as 
in key to symbols. Full line corresponds to w ~ t 1 ^ 3 . Inset: 
scaling plot of data on main diagram, using a = 1/2, z = 3/2; 
U = L/325. 



Figure 6: (Color online) Double-logarithmic plot of inter- 
face width against time, corresponding to TASEP with open 
boundary conditions inside the low-density (LD) phase, a/ = 
1/4, 0b = 1/2. Initial density p = 1/4, lattice sizes L as 
in key to symbols. Full line corresponds tom~ t 5 ' 12 . Inset: 
scaling plot of data on main diagram, using a = 1/2, z = 6/5; 
L, = L/325. 



multicritical point ai — (3e = 1/2, and found the KPZ 
exponents again there, with an accuracy similar to that 
obtained deep within the MC phase. The value z = 3/2 
has been found at this point by DMRG as well (l7l |. 

Elsewhere on the phase diagram, a low-density phase 
exists at ai < 1/2, ai < (3e (with p — a/), and a high- 
density phase at /3b < 1/2, Pe < cti (with p = 1 — 
Pe) UMMM- There is a critical coexistence line ai = 
(3e < 1/2, where a first-order transition occurs [l|. The 
low- and high-density phases are further subdivided 
However, such subdivisions will not concern us directly 
here, the relevant fact being that (outside the MC phase) 
the Bethe ansatz solution predicts a non-vanishing gap 
as L — > oo (i.e. z = 0) everywhere (this is found by 
DMRG as well [17]), except on the coexistence line where 
diffusive behavior with z = 2 is expected [18]. 

We first examined a point away from the coexistence 
line, namely ai — 1/4, (3e = 1/2 where the stationary 
density is thus p = 1/4. Setting the initial density at 
the stationary value, the time evolution of the associated 
interface widths is qualitatively very similar to the cases 
illustrated earlier, as shown in Fig. [6j however, scaling 
turns out to be very different. 

Attempting curve collapse on the data of Fig. [6] (see 
inset) gives the following estimates: a — 0.497(3), z — 
1.20(5). This would imply = 0.41(2) from scaling; a 
direct fit of L = 1300 data for 10 2 < t < 10 3 gives 
(3 = 0.35(1), which is a larger discrepancy compared with 



the scaling prediction than, e.g., for PBC or for the MC 
phase with open boundaries. 

These results for a and z are inconsistent with the scal- 
ing relation a + z = 2 from galilean invariance j2] ; how- 
ever, it should be recalled that translational invariance, 
the key ingredient for galilean invariance to hold, is in 
general broken by the system's boundaries here present. 
The fact that a + z = 2 is obeyed for open boundary 
conditions, in the MC phase a/, (3e > 1/2, can be ex- 
plained on the basis of simple kinematic- wave theory [l6j 
for the TASEP. Indeed, in this phase the kinematic waves 
produced by both boundaries do not penetrate the sys- 
tem [l3|, and one does not see the formation of a shock 
(density wave) in the bulk which would otherwise disrupt 
the translational symmetry. 

As mentioned above, both DMRG numerics and 
the Bethe ansatz solution [l8| predict z = 0, i.e., the 
correlation length is supposed to be finite here. We de- 
fer discussion, and a proposed solution, of the apparent 
contradiction between our own results and previous ones, 
to Subsection IIII Bl below, where a continuum mean-field 
treatment of the approach to stationarity is developed. 
For the moment, we note that an explanation can be pro- 
vided for the limiting- width exponent a, by going back to 
Eq. ([7]) which was there applied to the factorizable- weight 
case of PBC. As seen from the development immediately 
below Eq. J7]), the only changes brought about by a fi- 
nite correlation length (assumed to be <C L) amount to 
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Figure 7: (Color online) Double-logarithmic plot of inter- 
face width w and overall particle density p against time, 
corresponding to TASEP with open boundary conditions, at 
oti = Pe — 1/4 (on the coexistence line), lattice size L — 325, 
and two distinct initial densities, po = 1/4 and 1/2. 

an -L-independent correction to the term within round 
brackets, resulting from short-range density-density cor- 
relations. Thus, the L 1 / 2 dependence of [(w(L) 2 )] 
given in Eq. iJH), remains valid here. This is corrobo- 
rated by our numerical estimate a — 0.497(3), quoted 
above. 

Next, we looked at a point with olj = (3e = 1/4, on the 
coexistence line. In the stationary state, one expects a 
low-density phase with p_ = or, and a high-density one 
with p + = 1 — )3e, to coexist. We investigated the time 
evolution of both interface widths and overall densities, 
with two different initial conditions, namely po = 1/4, 
1/2 (because of particle-hole symmetry, starting with p = 
3/4 gives the same interface width as for p = 1/4, and 
a complementary particle density). It can be seen in 
Fig. [7] that with both initial conditions, the fixed point 
for the average density is p = 1/2. While, for po = 
1/4, the time evolution of the associated interface width 
still displays the simple, monotonic, character found in 
all setups previously studied here (apart from a small 
bump at early times), such a feature is lost for p = 
1/2. Though the width eventually settles at a unique 
saturation value, the corresponding relaxation time is one 
order of magnitude longer than elsewhere on the phase 
diagram or for PBC (see the curves for L — 325 in the 
respective figures). 

In order to unravel the corresponding spatial parti- 
cle distributions, we also looked at the time evolution of 
slope PDFs at this point, for the same initial densities. 

Results are depicted in Fig. [8j showing that for both 
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Figure 8: (Color online) Slope PDFs at assorted times, for 
interfaces corresponding to TASEP with open boundary con- 
ditions, at ai = Pe = 1/4 (on the coexistence line), lattice 
size L = 325, and two distinct initial densities, namely (a) 
po = 1/4, and (b) po = 1/2. Key to symbols: squares, t = 150 
(a), t = 50 (b); triangles, t = 500; hexagons, t = 5000; crosses, 
t = 20000. 



cases a double-peaked structure eventually evolves. The 
peak heights indicate a large degree of spatial segregation 
between the po = 1/4 and 3/4 phases. For example, in 
case (b) [initial density po = 1/2], at t = 20000 the peaks 
associated to the coexisting densities are ~ 8 times higher 
than the trough in the PDF at slope zero (p = 1/2). 

We pursued this point further, via direct examination 
of the evolution of averaged density profiles in the par- 
ticle system. Results for L = 325, with initial density 
po = 1/2 and initial distribution as uniform as possible 
(i.e., alternating empty and occupied sites), are shown in 
Fig.H 

For t < 500, the low injection/ejection rates cause the 
the initial plateau of uniform density to be symmetrically 
eaten into, as illustrated by the t = 50 profile. After 
the two density waves meet, a shock (kinematic wave) is 
formed, which for this case of ai — (3e < 1/2, is station- 
ary on average [16], i.e., it jiggles around and is in effect 
bounced off the system's boundaries. For the ensemble- 
averaged densities at fixed times > 500, the consequence 
of this is that the promediated profiles grow ever more 
featureless (because the spread between locations of the 
shock, at the same time but for different noise realiza- 
tions, increases with time owing to increasing sample-to- 
sample decorrelation) . At t — 10000, one sees in Fig. [9] 
a nearly-constant slope, i.e., the local average density in- 
creases roughly linearly with position along the system. 
Looking back at Fig. [7j however, it is apparent that a 
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Figure 9: (Color online) Ensemble-averaged particle density 
profiles at assorted times, for TASEP with open boundary 
conditions, at ai = (3e = 1/4 (on the coexistence line), lattice 
size L — 325, and initial density po = 1/2. Averages taken 
over 10 4 independent samples (see text). 



large degree of spatial segregation remains at stationar- 
ity between the p — 1/4 and 3/4 phases, within each 
realization. 

Finally, still on the coexistence line, we investigated the 
scaling properties of interface widths. For this, we chose 
the initial condition p — 1/4 which, as seen in Fig. [3 
results in a relatively smooth and monotonic evolution 
pattern. Our data (skipping the very early stage t < 10 3 , 
for which the small bumps referred to above make their 
appearance) are displayed in Fig. [TUl 

Attempting curve collapse on the data of Fig. [10] 
(see inset) gives the following estimates: a = 0.99(1), 
z = 2.10(5). This would imply (3 = 0.47(2) from scal- 
ing, as shown by the straight line on the main plot; a 
direct fit of L = 1300 data for 7 x 10 3 < t < 7 x 10 4 
gives (3 — 0.39(1), which again is a larger discrepancy 
compared with the scaling prediction than in cases pre- 
viously examined here. 

The estimate of z from curve collapse is roughly in 
line with the Bethe ansatz prediction f.18] z = 2, though 
it seems difficult to stretch the error bars for our data 
to include this latter value. As regards a, our estimate 
suggests that a = 1 is possibly an exact result for this 
case. 

An argument can be given as follows. Going back to 
the calculation outlined, for PBC, in Eq. ([7]), and recall- 
ing the spatial phase segregation shown in the late-time 
data of Fig. [HI one sees that the dominant feature of the 
interface (particle system) is its division at a strongly 
fluctuating interface into two main segments with sym- 



Figure 10: (Color online) Double-logarithmic plot of inter- 
face width against time, corresponding to TASEP with open 
boundary conditions on the coexistence line, at oti = [3e = 
1/4. Initial density po = 1/4, lattice sizes L as in key to sym- 
bols. Full line corresponds to w ~ t 0Ar . Inset: scaling plot of 
data on main diagram, using a = 1, z = 2.1; L» = L/325. 



metric slopes (low and high densities), of lengths oc L. 
In order words, local densities will be correlated along 
distances of order L. This is enough to guarantee that 
the mean square interface width will depend on L 2 . This 
argument is supported by a calculation of the size de- 
pendence of the stationary mean square width. On the 
coexistence line the diffusing shock gives a linear average 
stationary density profile (refer to the late-time curves in 
Fig. [9]), and this makes the large time limit of Eq. J7]) 
proportional to L 2 (2ai — l) 2 , consistent with a = 1. 



B. Continuum mean- field approach 

We investigate how the effects of number conservation 
can influence the approach to stationarity. We wish to 
find out whether, and under what conditions, transient 
(e.g., ballistic) phenomena can occur which would mask 
the underlying long-time relaxational dynamics. Deter- 
mination of the exponent z via interface width scaling, 
as done in Sees. ITTI and IIII A] relies heavily upon collaps- 
ing the 'shoulders' which mark the final approach to sta- 
tionarity, thus this method always picks out the longest 
characteristic time. Therefore one may ask, for the low- 
density phase with open boundary conditions, whether 
the L— independent relaxation time implied by the pre- 
diction z = is hidden underneath a longer process; and 
whether in this scenario the latter time, which is captured 
by interface width scaling, is associated with features not 
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so far emphasised within Bethe ansatz investigations, as 
in the connection between ballistic motion and imaginary 
parts of energies. Such things might explain the seeming 
mismatch between our results and those of Ref. [IE 

In the TASEP (and similar non-equilibrium flow pro- 
cesses) the bulk hopping process conserves total particle 
number - hence the continuity-equation nature of the evo- 
lution equation, and the slow nature of long-wavelength 
fluctuations of the density. With open boundaries, parti- 
cles can appear and leave at the chain ends, thus causing 
the total number of particles inside to change. 

A group of particles crossing at a boundary immedi- 
ately affects the mean density and the local density near 
the boundary. Such a group may be transferred into the 
interior as a kinematic wave/moving domain wall pro- 
vided the kinematic wave velocity v is non-zero and of 
the right sign. A time of order L/v will be needed to 
affect the local density throughout a system of size L (as 
for example is typically required to achieve a steady state 
configuration). 

This argument suggests that measurements of quanti- 
ties (like ones conserved in the interior) whose changes 
are boundary-induced will show characteristic times of 
order L z , with z' at least unity (still slower limiting pro- 
cesses may also be involved, owing to intrinsic features 
of the dynamics). 

In order to provide a quantitative counterpart to these 
ideas, we have used a continuum mean-field approach 
(see, e.g., Ref. 15) in which the system is described by 
the noiseless Burgers /KPZ equations, linearized by the 
Cole-Hopf transformation fig . [2p| . This picks up the 
kinematic wave effects in the evolution of the local den- 
sity p£ and its mean p. It captures the ballistic transport 
between boundary and interior, described qualitatively 
above, which gives z' = 1, In this description the re- 
duced density ae(t) — Pt(t) — 1/2 evolves like: 

a t (t) = ^ ln{coshif (I - £ ) e^ 2 ' + 

+Y^(A k e M + B k e- kt )e- k2t } , (10) 
k 

towards a steady state: 

at = K t&nhK(£-£ ) , (11) 

where K and £q are decided by boundary rates (a/, /3e), 
K being real, except in the MC phase where it is imag- 
inary. The sum takes care of the difference between 
initial and steady states, and (since this difference is de- 
caying) it involves complex k with Re(fc 2 ) > K 2 ; this 
k has to be consistent with boundary conditions. Af., 
Bk are therefore determined by initial conditions. The 
difference typically becomes a kinematic wave or do- 
main waii [a EE Hi HI HE HI HE Hi], if it is kmk- 

like, its velocity v depends as follows on coarse-grained 
densities p < , p> and currents J<, J> on either side of 
the kink: v — (J< — J > )/{p < — /»>). These quantities 
are set by the steady-state and initial (uniform) density 



profiles; for smooth small-amplitude kinematic waves 
v = dJ/dp = 1 — 2p, and then the local coarse-grained 
density p is at late times set by the steady state. 

For our investigations at (a>i, (3e) — (1/4, 1/2) the ini- 
tial profile pi = 1/4 has to go into the steady-state profile, 
which corresponds to if = 1/4, £ — L, and differs from 
the initial one only in having an upturn (to pi, = 3/8) at 
the right boundary. 

The kinematic wave velocity for p — 1/4 is 1/2, i.e., 
positive, so there is ballistic transfer, actually of vacan- 
cies at speed 1/2 from the right boundary, which after 
time t ~ L/v = 2L produces the steady-state profile with 
its upturn at the right boundary. This is corroborated 
by our numerical results depicted in Fig. [6] For all cases 
L = 325, 650, and 1300, the interface width at t = 2L has 
reached more than 98% of its asymptotic value. Thus, 
while attempts at producing overall curve collapse indi- 
cate z' ~ 1.20 as quoted above, by focusing exclusively 
on the scaling of the 'shoulders' one gets a result much 
closer to the mean-field value. 

We also applied this picture to other cases investigated 
in Subsection IIII A} to check for consistency. Results are 
as follows. 

For the MC phase at (a/, f3 E ) = (3/4,3/4), initially 
Pi = 1/2, so v = and the steady state pi differs from 
the average by an upturn to 2/3 at the left boundary, and 
a downturn to 1/3 at the right one. So ballistic effects 
are not very important; furthermore the limiting process 
is the slower "KPZ diffusion" having z = 3/2. 

On the coexistence line at (oli,[3e) = (1/4,1/4), two 
initial densities were used: (a) pi — 1/2, (b) pi = 1/4. 
For (a), the essential aspects are all as given in Fig.[9l i.e., 
a ballistic early evolution from the domain walls coming 
in from either end, followed by the limiting slower dif- 
fusive process (z — 2). These together give rise to the 
non-monotonic form in Fig. \7\ For (b) , the ballistic pro- 
cess is effective for a long time (of order L) during which 
the average p builds up, but again the limiting process is 
diffusion (z = 2). 

So, the mean-field continuum picture is consistent with 
all our numerical results for open boundary conditions, 
especially the one where we differ from Ref. EE In this 
latter case, the characteristic time arising from propaga- 
tion of the kinematic wave is longer than the intrinsic 
(L— independent) one, thus resulting in the effective ex- 
ponent z' = 1 (for all other cases, an exponent z > 1 as- 
sociated to intrinsic dynamics dominates anyway). Note 
that, for this scenario to work, one needs the kinematic 
wave to have non-zero effective velocity v, and (for a char- 
acteristic relaxation time, proportional to L, to show up) 
one also needs the kinematic wave to have to traverse a 
length of order L - whether that is the case depends on 
the maximum distance from the nearest effective bound- 
ary to the place(s) where the profile has to be adjusted. 

In summary, the difference between our results and 
those of Ref. Ha is real, and interpretable as above. The 
effects and interpretation may be important for other 
quantities conserved in the bulk which build up only by 
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propagation from the boundaries. 



IV. QUENCHED DISORDER 

Quenched disorder in the TASEP has been studied by 
many authors (27L [28l . I29I l30l. l3ll l32T| . Here, we investigate 
bond disorder [29], i.e., while all particles are identical, 
the site-to-site hopping rates are randomly distributed. 

Hereafter, we restrict ourselves to binary distributions 
P(p) for the internal nearest-neighbor hopping rate p: 



P(p) = tpS(p - p s ) + (1 - <p) 8(p - p w ) 



(12) 



where p s > 1/2, p w < 1/2 are associated, respectively, to 
"strong" and "weak" bonds. 

One must recall that the bond disorder introduced 
above gives rise to correlated, or "columnar" disorder 
in the associated interface problem [H, [33|. Indeed, the 
(fixed) value of the hopping rate at a given position x 
along the particle-model axis will determine, once and 
for all, the probability of the height h{x) being updated. 
This is in contrast with the usual picture of quenched 
disorder in the KPZ model El, [Hi, E3, EH, Ell , where 
it is assumed that the intensity of disorder is a random 
function of the instantaneous two-dimensional position 
(x, h(x)) of the interface element at x. 

In Figure [TT] we show interface width data for PBC, 
p = 1/2, and <p = 1/2, p s = 0.8, p w = 0.2 in Eq. (pjl. 
Note that the relaxation times are one order of magni- 
tude longer than is usual for pure systems (in the lat- 
ter case, one must make exception for the coexistence 
line with open boundary conditions). Though the over- 
all picture of a scaling regime still holds, with interface 
widths evolving in a simple, monotonic way, in general 
the quality of data collapse is lower than for pure systems. 
From our best fit, we estimate a = 1.05(5), z = 1.7(1), 
from which scaling gives (3 — 0.62(7). A direct fit of 
2 x 10 3 < x < 2 x 10 4 data for L = 1300 gives (3 = 0.56(1), 
only just within the error bars predicted by scaling. 

We have examined the stationary-state slope distribu- 
tions for the interface problem. Numerical results are 
displayed in Fig. [TJl together with their fit by a double- 
Gaussian form, <J>(x) = a G\(x) + (1 — a) (^(x), where Gi 
is a Gaussian centered at Xi with variance of. The best- 
fitting parameters give a roughly symmetric curve, with 
a = 0.53(2), x x = -0.26(1), x 2 = 0.31(1), a x = 0.21, 
(T2 = 0.18. Such a double-peaked structure indicates 
phase separation, as seen earlier for pure systems on 
the coexistence line (though here this is quantitatively 
milder, as the peak-to-trough ratio is ~ 1.3, to be com- 
pared to ~ 8 in the previous case). Phase separation 
is known to be a feature of the quenched-disordered 
TASEP 

We now outline a theoretical framework for the descrip- 
tion of the quenched-disordered problem in the spirit of 
earlier work by Tripathy and Barma 27|. The analytic 
understanding of slope distributions and, particulary, of 
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Figure 11: (Color online) Double- logarithmic plot of interface 
width against time, corresponding to TASEP with periodic 
boundary conditions and binary quenched disorder, with <p = 
1/2, p s = 0.8, p w = 0.2 [see Eq. CE])]. p = 1/2, lattice sizes 
L as in key to symbols. Full line corresponds to w ~ t ' 61 . 
Inset: scaling plot of data on main diagram, using a — 1.05, 
z = 1.7 (see text); L, = L/325. Each point is an average over 
10 000 independent realizations of quenched randomness. 
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Figure 12: (Color online) Slope PDF in stationary state, 
corresponding to TASEP with quenched disorder, ip — 1/2, 
p s = 0.8, p w = 0.2 [see Eq. {12])]. Lattice size L = 325 with 
PBC, p = 1/2. Full line is a double Gaussian fit to data 
(see text). 100 independent realizations of quenched random- 
ness were used, for each of which samples were collected from 
3 x 10 4 consecutive interface configurations. 
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currents and profiles, rests on the division of the sys- 
tem, for a given disorder configuration, into an alternat- 
ing succession of weak-bond and strong-bond segments 
(S w and S s ) each containing only bonds of one strength. 
The longest weak bond segments determine the current 
through the system. This is easily seen in a mean field 
account of the steady state of the particle system (with 
PBC). Here the constant current J yields the relation 
(profile map) 



Pt+i = 1 - 



J 



Pi Pt 



(13) 



where pt is the hopping rate from site £ to I + 1 , and pt 
is the mean occupation of site £. 

Throughout a weak bond segment S™ of length n 
the map involves the constant (reduced) current J /pt = 
J/p w = J w . The corresponding strong bond variable 
is J/p s = Js, Js < Jw since p s > p w - So, within 
any S™ (or S*) the profile map is that of an effec- 
tive pure system, which is well known to give density 
profiles of kink shape, corresponding to low current or 
high current: pt — \ — k taiihk(£ — £q) (monotonic 
increasing) or pt — \ = —KtanK(£ — £q) (monotonic 
decreasing), depending on whether the reduced current 
J/p is less than or greater than 1/4. k and K are re- 
lated to the reduced current by k = y/ (1/4) — (J/p); 

if=V(«)-(l/4) & 

In the high current case K has to be small [< 0(1/ L) 
in a segment of size L], to prevent the tangent from di- 
verging and taking pe outside of the permitted physi- 
cal range [0,1]. In the binary random system, it is not 
possible to have both weak and strong bond segments 
(having respectively K — K w , K s ) in the high current 
'state' since that would lead to monotonically decreas- 
ing p for all segments. That would violate the periodic 
boundary condition requirement. It cannot even apply 
with open boundary conditions in a large system, be- 
cause even with K s kept small by having J s close to 1/4, 
the larger J w would cause a large K w , resulting in non- 
physical values of pi. Arguments of this sort show that 
the strong bond segments S s are all in the low current 
phase, i.e., J s < 1/4 with density profiles in each 5 s in- 
creasing monotonically. Thus, for PBC the profiles in 
each S w have to decrease. That can come about from 
having J w > 1/4, which leads to K w = \/ J w — (1/4). 
If the longest weak bond segment has length no, to 
prevent unphysical pi's resulting from a divergence of 
ta,nK w (£ — £o) somewhere within that segment we must 
have K w < tt/uq. This gives: 



1 J , 1 / 7T 

- < — < - + — 

4 p w ~ 4 \n 



(14) 



It is straightforward to show that the characteristic 
length of weak bond segment is n ~ In (1/(1 — if)), 
and the largest weak segment length is n Q ~ 
ln_L/ln(l/(l — </?)), so the above condition on J w is very 
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Figure 13: (Color online) Upper part: time-averaged particle 
density against position on section of L = 325 system with 
PBC, corresponding to TASEP with quenched disorder, for a 
fixed bond configuration (seen in lower part of figure). Overall 
density p = 1/2, <p = 1/2, p s = 0.9, p w = 0.1 [see Eq. (p~2]) ] . 
Average is over 4.65 x 10 s successive time steps, in stationary 
state. Note kink-shaped profile for 35 < x < 45, in direct cor- 
respondence with the longest concentration of strong bonds; 
also, the high-current profile for 65 < x < 75, coinciding with 
the longest concentration of weak bonds. 



restrictive, making the current typically J ~ (p w /4) + 
0((l/lni) 2 ). 

This is in agreement with simulation results J ~ 0.055 
and 0.029, respectively for p w — 0.2 and p w = 0.1 (both 
with p — 1/2 on a lattice with L = 325). 

A decreasing profile in each S w is also possible with J w 
slightly less than 1/4 since, in addition to the kink-shaped 
low current profile, a profile with pt decreasing between 
1 and (1/2) + k can result from the low current map. 
Again the K w that allows that is limited, leading for this 
case to J ~ (Pmj/4) — O ((1/lnL) 2 ). This was seen, for 
example, in simulations for p = 0.8, where p w = 0.1 
gives J ~ 0.024. This same result can also arise from an- 
other low current profile (having pt decreasing between 
(1/2) — k and zero). All these weak segment profiles are 
seen in simulation results, together with the characteris- 
tic kink-shaped profiles of the strong-bond segments (e.g., 
for p w = 0.1, p = 0.5 see Fig. SH having J = 0.0297(3), 
i.e., J /p w just greater than 1/4). 

Local density or slope distributions P(v), with a = 
dh/dx, are available from the profiles pt, just as for 
the pure case, by using the known density-slope relation, 
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Eq. 10} in: 



y — 



segments 1 



dat/dl\ ^ \dai 

' ' segments 



(15) 

Contributions from high current segments give (using 
at = -KtsmKt): 



C = 



1 



2[a 2 +AK 2 ] 1 , (16) 



\dai/d£\ ae=a 
while low current segments have contributions 



_ f 2 [4fc 2 - d 2 ] for | a\ < 2k 
~ \2 [cr 2 -4fc 2 ] -1 for I cr| > 2k 



(17) 



In Eqs. lfl6]l . lfT7]). these contributions are superimposed 
with (,5-dependent weights related to their frequency of 
occurrence. Stochastic effects, missed from this account, 
cause fluctuations in the spatial pinning of the profiles 
(as for the pure case), leading mainly to a smoothing of 
the divergences near a = ±2k in the low current segment 
contributions. Quantitative comparison of these predic- 
tions, e.g., to the data displayed in Fig. [12], is not entirely 
appropriate, on account of the coarse-graining introduced 
by taking slope samples along segments of size m bonds 
(typically m = 60, as mentioned above). Nevertheless, 
one can see that some general aspects, such as the ap- 
pearance of peaks roughly equidistant from a = 0, are 
effectively mirrored in the numerical results. 



V. DISCUSSION AND CONCLUSIONS 

We start our discussion by recalling that the thermo- 
dynamic limits of interface models, such as KPZ, and 
of exclusion processes, differ in a subtle way. Namely, 
as can be seen from Eqs. (O and ([6]), if one takes the 
L — > 00 limit before allowing t — > 00, one is left with 
a perpetual coarsening transient in which the interface 
gets ever rougher. On the other hand, the TASEP dis- 
plays a well-defined stationary state even on an infinite 
system. Therefore, the correspondence of the limiting- 
width regime of the former type of problem to the sta- 
tionary state of the latter can be only take effect within 
a finite-size scaling context. 

Once one is mindful of this distinction, however, the 
similarities and differences between two finite systems 
linked by the correspondence recalled in Sec. [J are ex- 
pected to be bona fide features, which reflect objective 
connections between the physics of non-equilibrium flow 
processes, and that of moving interfaces in random me- 
dia. 

In Sec. HI] we first confirmed that the known KPZ ex- 
ponents can be numerically extracted, with good accu- 
racy, from the scaling of interface widths derived from 
the underlying TASEP with PBC. For a given range of 
system sizes, direct evaluation of (3 by examination of 



the transient regime of width growth against time seems 
to be the least accurate procedure, which in this case 
gave (3 — 0.31(1). From scaling, with a = 0.500(5), 
z = 1.52(3), one gets — 0.33(1), in much better agree- 
ment with the exact (3 = 1/3. Measuring this exponent 
directly from the transient regime tends to result in un- 
derestimation; as seen in Sees, [ill] and [IV] above, such a 
trend is present in all our subsequent results, both for 
open systems and for quenched disorder. Also for PBC, 
we showed that the dependence of limiting KPZ interface 
widths against TASEP particle density can be accounted 
for by a treatment, which makes explicit use of the weight 
factorization that occurs for TASEP with PBC Ar- 
guments based on weight factorization provide an expla- 
nation for the shape of slope distributions in the interface 
problem as well. 

In Subsec. MI A[ we first established that open- 
boundary systems in the maximal-current phase, charac- 
terized by ai, [3e > 1/2, exhibit the same set of KPZ 
exponents as their PBC counterparts. For the expo- 
nent z. this is in agreement with the Bethe ansatz so- 
lution [18| ■ Examination of interface widths correspond- 
ing to systems in the low-density phase, with ai = 1/4, 
(3e = 1/2, shows that curve collapse can be found to a 
rather good extent, giving the following exponent esti- 
mates: a = 0.497(3), z = 1.20(5), j3 = 0.41(2). This 
apparent disagreement with the Bethe ansatz predic- 
tion [l8| of z = 0, which implies a finite correlation 
length, is addressed in Subsec. IIIIBI (see two paragraphs 
on). We have been able to provide a prediction of the 
(possibly exact) interface width exponent a = 1/2, based 
on considerations which make explicit use of a finite cor- 
relation length for the TASEP. 

For systems with ai = (3e < 1/2, i.e., on the coexis- 
tence line of the open-boundary TASEP, interface width 
scaling gave a = 0.99(1), z = 2.10(5), (3 = 0.47(2). We 
provided a prediction of a — 1 based on properties of the 
corresponding TASEP, in this case the fact that phase 
separation governs the dominant features of the interface 
configuration at stationarity. Direct calculation of den- 
sity profiles in the particle system shows the time evo- 
lution of a shock (kinematic wave). At late times, the 
ensemble-averaging of local densities tends to mask the 
evidence of phase segregation, which can, however, be 
retrieved by examination of the corresponding interface 
slope PDFs. 

In Subsec. IIIIB[ we outlined a mean-field contin- 
uum calculation, which sheds additional light on the ap- 
proach to stationarity in open-boundary systems. We 
showed that, under suitable conditions such as those 
at (ui,[3e) = (1/4,1/2) with uniform initial density 
pi = 1/4, system- wide propagation of a kinematic wave 
translates into a characteristic time « I 2 (z' = 1 in 
mean field). This goes towards explaining the appar- 
ent inconsistency between our result from interface- width 
evolution, z = 1.20(5), and that from the Bethe ansatz 
solution which gives z = 0. Indeed, in that case the 
L— independent relaxation time implied by z = is hid- 
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den underneath a slower part of ballistic origin, and it 
is the scaling of the latter which is captured by the in- 
terface width collapse, but not by considerations solely 
of the real part of Bethe ansatz excitation energies. Fur- 
thermore, in all other cases where our own numerical 
results are consistent with those of Ref. the mean- 
field prediction concurs with both. An extension of this 
work has now been carried out using the more precise 
domain wall method of Refs .[2p . [40l on a kink-like initial 
state in the massive phase [41j. This shows clearly the 
ballistic element and the much faster (size-independent) 
amplitude decay, providing an independent confirmation 
of the scenario. 

In Sec. [IV] we investigated quenched bond disorder in 
TASEP with PBC. From the scaling analysis of the cor- 
responding interface widths (which in this case are sub- 
jected to correlated, or "columnar" randomness [29l. l33|) 
we estimate a = 1.05(5), z = 1.7(1), (3 = 0.62(7). 
For comparison, values quoted for (3 in standard, two- 
dimensional quenched disorder in D = 1 + 1 KPZ systems 
are close to 0.63 [HEH]. It has been argued (HEa] that 
this type of KPZ model is in the universality class of 
directed percolation, thus one should have (3 = (3p,p = 
0.633. Furthermore, in Ref. 1371 it is shown that by vary- 
ing the intensity of the various terms in the quenched 
counterpart of Eq. ([3]), one can make KPZ-like sys- 
tems go through distinct regimes, namely: pinned, with 
a P = 0.63(3), (3 P = 0.67(5), z P = 1.06(8); moving, with 
OL M = 0.75(4), [3 M = 0.74(6), z M = 1.01(10); and an- 
nealed (i.e., fast-moving interface), with cca = 0.50(4), 
(3a = 0.30(4), za = 1.67(26). It can be seen that our 
own set of estimates does not fully fit into any of these, 
as could reasonably be expected from the extreme cor- 
relation between quenched defects which is present here, 
and not in those early examples. 

We tested universality properties within the columnar 
disorder class of models. This was done by replacing the 
binary distribution, Eq. lfl~2|). with a continuous, uniform 
one: P(p) = (1 — c) _1 for c < p < 1. We used c = 0.1, 
thereby avoiding the problems associated with allowing 
p = 0, see Ref. [23, while still having a rather broad dis- 
tribution. Our data scale similarly to those displayed in 
FigureQTjfor Eq. H2). We get a = 1.05(5), z = 1.45(10) 
(hence (3 — 0.72(7) from scaling), in reasonable agree- 
ment with the binary disorder case, though error bars 
for z just fail to overlap. 

We also investigated slope distributions for the 
quenched disorder problem. These provide clear evidence 
of phase separation, a phenomenon known to take place 



in such circumstances [27]. 

A direct analytic (mean-field) approach to steady-state 
properties of TASEP with quenched disorder produced 
closed-form expressions for the piecewise shape of aver- 
aged profiles densities, as well as rather restrictive bounds 
on currents. All these have been verified in our numeri- 
cal simulations. The analytic approach is similar to that 
of Ref. [27], where it was already shown that a mean field 
description applies, and to part of Ref. 132. In place of the 
maximum current principle used in Ref. [271 we have ob- 
tained analytic consequences of the mean field mapping 
within segments and combined them with segment prob- 
abilities to obtain new results, particularly for profiles 
and limiting currents. 

We note that, for weak randomness, characterized by 
a small value of a disorder parameter e [this could be, 
e.g., p s —p w in Eq. JUJ], and steady state conditions, the 
noiseless (mean-field) constant current condition gives an 
equation for a = dh/dx in the form: 

ec( * )= £ +fT2 ' (i8) 

where ((x) is a quenched variable corresponding to ran- 
dom bond disorder. Eq. fl8j) has the same structure as 
that found for equations governing the evolution of lo- 
cal Lyapunov exponents for Heisenberg-Mattis spin glass 
chains [H, S3|- To see the correspondence, refer, e.g., to 
Equation (14) of Ref. 0, substituting e for (low) magnon 
frequency u>. So one can, in principle, adopt the same 
type of Fokker-Planck procedures to find distributions of 
a in Burgers-like equations, and hence of slope distribu- 
tions in KPZ systems. 
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